
library(Matching)
library(tidyverse)

################################################################################

rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. Dryland: peat depth = 0

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW==0)

# did outcome vars

virginclip_pop$did2015 <- virginclip_pop$forestha_2015 - virginclip_pop$forestha_2011
virginclip_pop$did2014 <- virginclip_pop$forestha_2014 - virginclip_pop$forestha_2011
virginclip_pop$did2013 <- virginclip_pop$forestha_2013 - virginclip_pop$forestha_2011
virginclip_pop$did2012 <- virginclip_pop$forestha_2012 - virginclip_pop$forestha_2011

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

# time indices 

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2015
t3 = 2014
t4 = 2013
t5 = 2012

# ddd outcome vars

virginclip_pop$ddd2015 <- virginclip_pop$did2015 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2014 <- virginclip_pop$did2014 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2013 <- virginclip_pop$did2013 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2012 <- virginclip_pop$did2012 - ((t5-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean +  pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

# attach propensity score to data
virginclip_pop$pscore<- ps_virginclip_pop$fitted

##### DID REGRESSIONS #####

PSM_DID_virginclip_pop2015 <- Match( Y = virginclip_pop$did2015, Tr = virginclip_pop$moratorium,
                              X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                              replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2015)



PSM_DID_virginclip_pop2014 <- Match( Y = virginclip_pop$did2014, Tr = virginclip_pop$moratorium,
                             X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                             replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2014)

PSM_DID_virginclip_pop2013 <- Match( Y = virginclip_pop$did2013, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2013)

PSM_DID_virginclip_pop2012 <- Match( Y = virginclip_pop$did2012, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2012)

##### DDD REGRESSIONS #####

PSM_DDD_virginclip_pop2015 <- Match( Y = virginclip_pop$ddd2015, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2015)

PSM_DDD_virginclip_pop2014 <- Match( Y = virginclip_pop$ddd2014, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2014)


PSM_DDD_virginclip_pop2013 <- Match( Y = virginclip_pop$ddd2013, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2013)

PSM_DDD_virginclip_pop2012 <- Match( Y = virginclip_pop$ddd2012, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2012)

# remove data and save results

rm(virginclip_pop, ps_virginclip_pop)
save.image("~/results/dryland-np-did-ddd-1215.RData")


### Peatland ###

rm(list=ls())
load("np-ddd/virginclip_pop.RData")

################################################################################

# 2. Peatland: peat depth > 0 

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0)

# did outcome vars

virginclip_pop$did2015 <- virginclip_pop$forestha_2015 - virginclip_pop$forestha_2011
virginclip_pop$did2014 <- virginclip_pop$forestha_2014 - virginclip_pop$forestha_2011
virginclip_pop$did2013 <- virginclip_pop$forestha_2013 - virginclip_pop$forestha_2011
virginclip_pop$did2012 <- virginclip_pop$forestha_2012 - virginclip_pop$forestha_2011

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

# time indices

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2015
t3 = 2014
t4 = 2013
t5 = 2012

# ddd outcome vars

virginclip_pop$ddd2015 <- virginclip_pop$did2015 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2014 <- virginclip_pop$did2014 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2013 <- virginclip_pop$did2013 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2012 <- virginclip_pop$did2012 - ((t5-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean +  pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

# attach propensity score to data
virginclip_pop$pscore<- ps_virginclip_pop$fitted

##### DID REGRESSIONS #####

PSM_DID_virginclip_pop2015 <- Match( Y = virginclip_pop$did2015, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2015)



PSM_DID_virginclip_pop2014 <- Match( Y = virginclip_pop$did2014, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2014)

PSM_DID_virginclip_pop2013 <- Match( Y = virginclip_pop$did2013, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2013)

PSM_DID_virginclip_pop2012 <- Match( Y = virginclip_pop$did2012, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2012)


##### DDD REGRESSIONS #####

PSM_DDD_virginclip_pop2015 <- Match( Y = virginclip_pop$ddd2015, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2015)

PSM_DDD_virginclip_pop2014 <- Match( Y = virginclip_pop$ddd2014, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2014)


PSM_DDD_virginclip_pop2013 <- Match( Y = virginclip_pop$ddd2013, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2013)

PSM_DDD_virginclip_pop2012 <- Match( Y = virginclip_pop$ddd2012, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2012)

# remove data and save results

rm(virginclip_pop, ps_virginclip_pop)
save.image("~/results/peatland-np-did-ddd-1215.RData")